home *** CD-ROM | disk | FTP | other *** search
/ Aminet 35 / Aminet 35 (2000)(Schatztruhe)[!][Feb 2000].iso / Aminet / gfx / misc / gnuplot-src.lha / gnuplot-3.7.1src / gnuplot-3.7.1.lha / gnuplot-3.7.1 / plot3d.c < prev    next >
Encoding:
C/C++ Source or Header  |  1998-12-10  |  46.8 KB  |  1,577 lines
  1. #ifndef lintnew_title, xp->tSid = "$Id: plot3d.c,v 1.16 1998/12/10 18:30:52 lhecking Exp $";
  2. #endif
  3.  
  4. /* GNUPLOT - plot3d.c */
  5.  
  6. /*[
  7.  * Copyright 1986 - 1993, 1998   Thomas Williams, Colin Kelley
  8.  *
  9.  * Permission to use, copy, and distribute this software and its
  10.  * documentation for any purpose with or without fee is hereby granted,
  11.  * provided that the above copyright notice appear in all copies and
  12.  * that both that copyright notice and this permission notice appear
  13.  * in supporting documentation.
  14.  *
  15.  * Permission to modify the software is granted, but not the right to
  16.  * distribute the complete modified source code.  Modifications are to
  17.  * be distributed as patches to the released version.  Permission to
  18.  * distribute binaries produced by compiling modified sources is granted,
  19.  * provided you
  20.  *   1. distribute the corresponding source modifications from the
  21.  *    released version in the form of a patch file along with the binaries,
  22.  *   2. add special version identification to distinguish your version
  23.  *    in addition to the base release version number,
  24.  *   3. provide your name and address as the primary contact for the
  25.  *    support of your modified version, and
  26.  *   4. retain our contact information in regard to use of the base
  27.  *    software.
  28.  * Permission to distribute the released version of the source code along
  29.  * with corresponding source modifications in the form of a patch file is
  30.  * granted with same provisions 2 through 4 for binary distributions.
  31.  *
  32.  * This software is provided "as is" without express or implied warranty
  33.  * to the extent permitted by applicable law.
  34. ]*/
  35.  
  36. #include "plot.h"
  37. #include "setshow.h"
  38. #include "binary.h"
  39.  
  40. #ifndef _Windows
  41. # include "help.h"
  42. #endif
  43.  
  44. #ifndef STDOUT
  45. #define STDOUT 1
  46. #endif
  47.  
  48. /* static prototypes */
  49.  
  50. static void get_3ddata __PROTO((struct surface_points * this_plot));
  51. static void print_3dtable __PROTO((int pcount));
  52. static void eval_3dplots __PROTO((void));
  53. static void grid_nongrid_data __PROTO((struct surface_points * this_plot));
  54. static void parametric_3dfixup __PROTO((struct surface_points * start_plot, int *plot_num));
  55.  
  56. /* the curves/surfaces of the plot */
  57. struct surface_points *first_3dplot = NULL;
  58. static struct udft_entry plot_func;
  59.  
  60.  
  61. extern struct udft_entry *dummy_func;
  62. extern int datatype[];
  63. extern char timefmt[];
  64. extern TBOOLEAN is_3d_plot;
  65. extern int plot_token;
  66.  
  67. /* in order to support multiple axes, and to
  68.  * simplify ranging in parametric plots, we use
  69.  * arrays to store some things. For 2d plots,
  70.  * elements are  y1 = 0 x1 = 1 y2 = 2 x2 = 3
  71.  * for 3d,  z = 0, x = 1, y = 2
  72.  * these are given symbolic names in plot.h
  73.  */
  74.  
  75. extern double min_array[AXIS_ARRAY_SIZE], max_array[AXIS_ARRAY_SIZE];
  76. extern int auto_array[AXIS_ARRAY_SIZE];
  77. extern TBOOLEAN log_array[AXIS_ARRAY_SIZE];
  78. extern double base_array[AXIS_ARRAY_SIZE];
  79. extern double log_base_array[AXIS_ARRAY_SIZE];
  80.  
  81. /* some file-wide variables to store which axis we are using */
  82. static int x_axis, y_axis, z_axis;
  83.  
  84.  
  85. /* Deleted from setshow.h and renamed */
  86. extern FILE *gpoutfile;
  87.  
  88. /* info from datafile module */
  89. extern int df_datum;
  90. extern int df_line_number;
  91. extern int df_no_use_specs;
  92. extern int df_eof;
  93. extern int df_timecol[];
  94. extern TBOOLEAN df_matrix;
  95.  
  96. #define Inc_c_token if (++c_token >= num_tokens)    \
  97.                         int_error ("Syntax error", c_token);
  98.  
  99. /* From plot2d.c */
  100. extern int reverse_range[];
  101.  
  102. /*
  103.  * IMHO, code is getting too cluttered with repeated chunks of
  104.  * code. Some macros to simplify, I hope.
  105.  *
  106.  * do { } while(0) is comp.lang.c recommendation for complex macros
  107.  * also means that break can be specified as an action, and it will
  108.  * 
  109.  */
  110.  
  111. /*  copy scalar data to arrays
  112.  * optimiser should optimise infinite away
  113.  * dont know we have to support ranges [10:-10] - lets reverse
  114.  * it for now, then fix it at the end.
  115.  */
  116. #define INIT_ARRAYS(axis, min, max, auto, is_log, base, log_base, infinite) \
  117. do{if ((auto_array[axis] = auto) == 0 && max<min) {\
  118.     min_array[axis] = max;\
  119.    max_array[axis] = min; /* we will fix later */ \
  120.  } else { \
  121.     min_array[axis] = (infinite && (auto&1)) ? VERYLARGE : min; \
  122.     max_array[axis] = (infinite && (auto&2)) ? -VERYLARGE : max; \
  123.  } \
  124.  log_array[axis] = is_log; base_array[axis] = base; log_base_array[axis] = log_base;\
  125. }while(0)
  126.  
  127. /* handle reversed ranges */
  128. #define CHECK_REVERSE(axis) \
  129. do{\
  130.  if (auto_array[axis] == 0 && max_array[axis] < min_array[axis]) {\
  131.   double temp = min_array[axis]; min_array[axis] = max_array[axis]; max_array[axis] = temp;\
  132.   reverse_range[axis] = 1; \
  133.  } else reverse_range[axis] = (range_flags[axis]&RANGE_REVERSE); \
  134. }while(0)
  135.  
  136.  
  137. /* get optional [min:max] */
  138. #define LOAD_RANGE(axis) \
  139. do {\
  140.  if (equals(c_token, "[")) { \
  141.   c_token++; \
  142.   auto_array[axis] = load_range(axis,&min_array[axis], &max_array[axis], auto_array[axis]);\
  143.   if (!equals(c_token, "]"))\
  144.    int_error("']' expected", c_token);\
  145.   c_token++;\
  146.  }\
  147. } while (0)
  148.  
  149.  
  150. /* store VALUE or log(VALUE) in STORE, set TYPE as appropriate
  151.  * Do OUT_ACTION or UNDEF_ACTION as appropriate
  152.  * adjust range provided type is INRANGE (ie dont adjust y if x is outrange
  153.  * VALUE must not be same as STORE
  154.  */
  155.  
  156. #define STORE_WITH_LOG_AND_FIXUP_RANGE(STORE, VALUE, TYPE, AXIS, OUT_ACTION, UNDEF_ACTION)\
  157. do { if (log_array[AXIS]) { if (VALUE<0.0) {TYPE = UNDEFINED; UNDEF_ACTION; break;} \
  158.               else if (VALUE == 0.0){STORE = -VERYLARGE; TYPE = OUTRANGE; OUT_ACTION; break;} \
  159.               else { STORE = log(VALUE)/log_base_array[AXIS]; } \
  160.      } else STORE = VALUE; \
  161.      if (TYPE != INRANGE) break;  /* dont set y range if x is outrange, for example */ \
  162.      if ( VALUE<min_array[AXIS] ) { \
  163.       if (auto_array[AXIS] & 1) min_array[AXIS] = VALUE; else { TYPE = OUTRANGE; OUT_ACTION; break; }  \
  164.      } \
  165.      if ( VALUE>max_array[AXIS] ) { \
  166.       if (auto_array[AXIS] & 2) max_array[AXIS] = VALUE; else { TYPE = OUTRANGE; OUT_ACTION; }   \
  167.      } \
  168. } while(0)
  169.  
  170. /* use this instead empty macro arguments to work around NeXT cpp bug */
  171. /* if this fails on any system, we might use ((void)0) */
  172. #define NOOP            /* */
  173.  
  174. /* check range and take logs of min and max if logscale
  175.  * this also restores min and max for ranges like [10:-10]
  176.  */
  177.  
  178. #ifdef HAVE_STRINGIZE
  179. # define RANGE_MSG(x) #x " range is less than threshold : see `set zero`"
  180. #else
  181. # define RANGE_MSG(x) "x range is less than threshold : see `set zero`"
  182. #endif
  183.  
  184. #define FIXUP_RANGE_FOR_LOG(AXIS, WHICH) \
  185. do { if (reverse_range[AXIS]) { \
  186.       double temp = min_array[AXIS]; \
  187.       min_array[AXIS] = max_array[AXIS]; \
  188.       max_array[AXIS] = temp; \
  189.      }\
  190.      if (log_array[AXIS]) { \
  191.       if (min_array[AXIS] <= 0.0 || max_array[AXIS] <= 0.0) \
  192.        int_error(RANGE_MSG(WHICH), NO_CARET); \
  193.       min_array[AXIS] = log(min_array[AXIS])/log_base_array[AXIS]; \
  194.       max_array[AXIS] = log(max_array[AXIS])/log_base_array[AXIS];  \
  195.     } \
  196. } while(0)
  197.  
  198.  
  199.  
  200. /* support for dynamic size of input line */
  201.  
  202. void plot3drequest()
  203. /*
  204.  * in the parametric case we would say splot [u= -Pi:Pi] [v= 0:2*Pi] [-1:1]
  205.  * [-1:1] [-1:1] sin(v)*cos(u),sin(v)*cos(u),sin(u) in the non-parametric
  206.  * case we would say only splot [x= -2:2] [y= -5:5] sin(x)*cos(y)
  207.  * 
  208.  */
  209. {
  210.     TBOOLEAN changed;
  211.     int dummy_token0 = -1, dummy_token1 = -1;
  212.  
  213.     is_3d_plot = TRUE;
  214.  
  215.     if (parametric && strcmp(dummy_var[0], "t") == 0) {
  216.     strcpy(dummy_var[0], "u");
  217.     strcpy(dummy_var[1], "v");
  218.     }
  219.     autoscale_lx = autoscale_x;
  220.     autoscale_ly = autoscale_y;
  221.     autoscale_lz = autoscale_z;
  222.  
  223.     if (!term)            /* unknown */
  224.     int_error("use 'set term' to set terminal type first", c_token);
  225.  
  226.     if (equals(c_token, "[")) {
  227.     c_token++;
  228.     if (isletter(c_token)) {
  229.         if (equals(c_token + 1, "=")) {
  230.         dummy_token0 = c_token;
  231.         c_token += 2;
  232.         } else {
  233.         /* oops; probably an expression with a variable. */
  234.         /* Parse it as an xmin expression. */
  235.         /* used to be: int_error("'=' expected",c_token); */
  236.         }
  237.     }
  238.     changed = parametric ? load_range(U_AXIS, &umin, &umax, autoscale_lu) : load_range(FIRST_X_AXIS, &xmin, &xmax, autoscale_lx);
  239.     if (!equals(c_token, "]"))
  240.         int_error("']' expected", c_token);
  241.     c_token++;
  242.     /* if (changed) */
  243.     if (parametric)
  244.         /* autoscale_lu = FALSE; */
  245.         autoscale_lu = changed;
  246.     else
  247.         /* autoscale_lx = FALSE; */
  248.         autoscale_lx = changed;
  249.     }
  250.     if (equals(c_token, "[")) {
  251.     c_token++;
  252.     if (isletter(c_token)) {
  253.         if (equals(c_token + 1, "=")) {
  254.         dummy_token1 = c_token;
  255.         c_token += 2;
  256.         } else {
  257.         /* oops; probably an expression with a variable. */
  258.         /* Parse it as an xmin expression. */
  259.         /* used to be: int_error("'=' expected",c_token); */
  260.         }
  261.     }
  262.     changed = parametric ? load_range(V_AXIS, &vmin, &vmax, autoscale_lv) : load_range(FIRST_Y_AXIS, &ymin, &ymax, autoscale_ly);
  263.     if (!equals(c_token, "]"))
  264.         int_error("']' expected", c_token);
  265.     c_token++;
  266.     /* if (changed) */
  267.     if (parametric)
  268.         /* autoscale_lv = FALSE; */
  269.         autoscale_lv = changed;
  270.     else
  271.         /* autoscale_ly = FALSE; */
  272.         autoscale_ly = changed;
  273.     }
  274.     if (parametric) {
  275.     /* set optional x (parametric) or z ranges */
  276.     if (equals(c_token, "[")) {
  277.         c_token++;
  278.         autoscale_lx = load_range(FIRST_X_AXIS, &xmin, &xmax, autoscale_lx);
  279.         if (!equals(c_token, "]"))
  280.         int_error("']' expected", c_token);
  281.         c_token++;
  282.     }
  283.     /* set optional y ranges */
  284.     if (equals(c_token, "[")) {
  285.         c_token++;
  286.         autoscale_ly = load_range(FIRST_Y_AXIS, &ymin, &ymax, autoscale_ly);
  287.         if (!equals(c_token, "]"))
  288.         int_error("']' expected", c_token);
  289.         c_token++;
  290.     }
  291.     }                /* parametric */
  292.     if (equals(c_token, "[")) {    /* set optional z ranges */
  293.     c_token++;
  294.     autoscale_lz = load_range(FIRST_Z_AXIS, &zmin, &zmax, autoscale_lz);
  295.     if (!equals(c_token, "]"))
  296.         int_error("']' expected", c_token);
  297.     c_token++;
  298.     }
  299.     CHECK_REVERSE(FIRST_X_AXIS);
  300.     CHECK_REVERSE(FIRST_Y_AXIS);
  301.     CHECK_REVERSE(FIRST_Z_AXIS);
  302.  
  303.     /* use the default dummy variable unless changed */
  304.     if (dummy_token0 >= 0)
  305.     copy_str(c_dummy_var[0], dummy_token0, MAX_ID_LEN);
  306.     else
  307.     (void) strcpy(c_dummy_var[0], dummy_var[0]);
  308.  
  309.     if (dummy_token1 >= 0)
  310.     copy_str(c_dummy_var[1], dummy_token1, MAX_ID_LEN);
  311.     else
  312.     (void) strcpy(c_dummy_var[1], dummy_var[1]);
  313.  
  314.     eval_3dplots();
  315. }
  316.  
  317.  
  318. static void grid_nongrid_data(this_plot)
  319. struct surface_points *this_plot;
  320. {
  321.     int i, j, k;
  322.     double x, y, z, w, dx, dy, xmin, xmax, ymin, ymax;
  323.     struct iso_curve *old_iso_crvs = this_plot->iso_crvs;
  324.     struct iso_curve *icrv, *oicrv, *oicrvs;
  325.  
  326.     /* Compute XY bounding box on the original data. */
  327.     xmin = xmax = old_iso_crvs->points[0].x;
  328.     ymin = ymax = old_iso_crvs->points[0].y;
  329.     for (icrv = old_iso_crvs; icrv != NULL; icrv = icrv->next) {
  330.     struct coordinate GPHUGE *points = icrv->points;
  331.  
  332.     for (i = 0; i < icrv->p_count; i++, points++) {
  333.         if (xmin > points->x)
  334.         xmin = points->x;
  335.         if (xmax < points->x)
  336.         xmax = points->x;
  337.         if (ymin > points->y)
  338.         ymin = points->y;
  339.         if (ymax < points->y)
  340.         ymax = points->y;
  341.     }
  342.     }
  343.  
  344.     dx = (xmax - xmin) / (dgrid3d_col_fineness - 1);
  345.     dy = (ymax - ymin) / (dgrid3d_row_fineness - 1);
  346.  
  347.     /* Create the new grid structure, and compute the low pass filtering from
  348.      * non grid to grid structure.
  349.      */
  350.     this_plot->iso_crvs = NULL;
  351.     this_plot->num_iso_read = dgrid3d_col_fineness;
  352.     this_plot->has_grid_topology = TRUE;
  353.     for (i = 0, x = xmin; i < dgrid3d_col_fineness; i++, x += dx) {
  354.     struct coordinate GPHUGE *points;
  355.  
  356.     icrv = iso_alloc(dgrid3d_row_fineness + 1);
  357.     icrv->p_count = dgrid3d_row_fineness;
  358.     icrv->next = this_plot->iso_crvs;
  359.     this_plot->iso_crvs = icrv;
  360.     points = icrv->points;
  361.  
  362.     for (j = 0, y = ymin; j < dgrid3d_row_fineness; j++, y += dy, points++) {
  363.         z = w = 0.0;
  364.  
  365. #ifndef BUGGY_DGRID_RANGING /* HBB 981209 */
  366.         /* as soon as ->type is changed to UNDEFINED, break out of
  367.          * two inner loops! */
  368.         points->type = INRANGE;
  369. #endif
  370.         for (oicrv = old_iso_crvs; oicrv != NULL; oicrv = oicrv->next) {
  371.         struct coordinate GPHUGE *opoints = oicrv->points;
  372.         for (k = 0; k < oicrv->p_count; k++, opoints++) {
  373.             double dist, dist_x = fabs(opoints->x - x), dist_y = fabs(opoints->y - y);
  374.  
  375.             switch (dgrid3d_norm_value) {
  376.             case 1:
  377.             dist = dist_x + dist_y;
  378.             break;
  379.             case 2:
  380.             dist = dist_x * dist_x + dist_y * dist_y;
  381.             break;
  382.             case 4:
  383.             dist = dist_x * dist_x + dist_y * dist_y;
  384.             dist *= dist;
  385.             break;
  386.             case 8:
  387.             dist = dist_x * dist_x + dist_y * dist_y;
  388.             dist *= dist;
  389.             dist *= dist;
  390.             break;
  391.             case 16:
  392.             dist = dist_x * dist_x + dist_y * dist_y;
  393.             dist *= dist;
  394.             dist *= dist;
  395.             dist *= dist;
  396.             break;
  397.             default:
  398.             dist = pow(dist_x, (double) dgrid3d_norm_value) +
  399.                 pow(dist_y, (double) dgrid3d_norm_value);
  400.             break;
  401.             }
  402.  
  403.             /* The weight of this point is inverse proportional
  404.              * to the distance.
  405.              */
  406.             if (dist == 0.0) {
  407. #ifndef BUGGY_DGRID_RANGING
  408.                 /* HBB 981209: revised flagging as undefined */
  409.             /* Supporting all those infinities on various
  410.              * platforms becomes tiresome, to say the least :-(
  411.              * Let's just return the first z where this happens,
  412.              * unchanged, and be done with this, period. */
  413.                 points->type = UNDEFINED;
  414.             z = opoints->z;
  415.             w = 1.0;
  416.             break; /* out of for (k...) loop */
  417. #else
  418. #if !defined(AMIGA_SC_6_1) && !defined(__PUREC__)
  419.             dist = VERYLARGE;
  420. #else /* !AMIGA_SC_6_1 && !__PUREC__ */
  421.             /* Multiplying VERYLARGE by opoints->z below
  422.              * might yield Inf (i.e. a number that can't
  423.              * be represented on the machine). This will
  424.              * result in points->z being set to NaN. It's
  425.              * better to have a pretty large number that is
  426.              * also on the safe side... The numbers that are
  427.              * read by gnuplot are float values anyway, so
  428.              * they can't be bigger than FLT_MAX. So setting
  429.              * dist to FLT_MAX^2 will make dist pretty large
  430.              * with respect to any value that has been read. */
  431.             dist = ((double) FLT_MAX) * ((double) FLT_MAX);
  432. #endif /* !AMIGA_SC_6_1 && !__PUREC__ */
  433. #endif /* BUGGY_DGRID_RANGING */
  434.             } else
  435.             dist = 1.0 / dist;
  436.  
  437.             z += opoints->z * dist;
  438.             w += dist;
  439.         }
  440. #ifndef BUGGY_DGRID_RANGING
  441.         if (points->type != INRANGE)
  442.           break; /* out of the second-inner loop as well ... */
  443. #endif
  444.         }
  445.  
  446. #ifndef BUGGY_DGRID_RANGING
  447.         /* Now that we've escaped the loops safely, we know that we
  448.          * do have a good value in z and w, so we can proceed just as
  449.          * if nothing had happened at all. Nice, isn't it? */
  450.         points->type = INRANGE;
  451.  
  452.         STORE_WITH_LOG_AND_FIXUP_RANGE(points->x, x, points->type, x_axis, NOOP, continue);
  453.         STORE_WITH_LOG_AND_FIXUP_RANGE(points->y, y, points->type, y_axis, NOOP, continue);
  454.         STORE_WITH_LOG_AND_FIXUP_RANGE(points->z, z / w, points->type, z_axis, NOOP, continue);
  455. #else
  456.         /* HBB 981026: original, short version of this code */
  457.         points->x = x;
  458.         points->y = y;
  459.         points->z = z / w;
  460.         points->type = INRANGE;
  461. #endif
  462.     }
  463.     }
  464.  
  465.     /* Delete the old non grid data. */
  466.     for (oicrvs = old_iso_crvs; oicrvs != NULL;) {
  467.     oicrv = oicrvs;
  468.     oicrvs = oicrvs->next;
  469.     iso_free(oicrv);
  470.     }
  471. }
  472.  
  473. static void get_3ddata(this_plot)
  474. struct surface_points *this_plot;
  475. /* this_plot->token is end of datafile spec, before title etc
  476.  * will be moved passed title etc after we return
  477.  */
  478. {
  479.     int xdatum = 0;
  480.     int ydatum = 0;
  481.     int i, j;
  482.     double v[3];
  483.     int pt_in_iso_crv = 0;
  484.     struct iso_curve *this_iso;
  485.  
  486.     if (mapping3d == MAP3D_CARTESIAN) {
  487.     if (df_no_use_specs == 2)
  488.         int_error("Need 1 or 3 columns for cartesian data", this_plot->token);
  489.     } else {
  490.     if (df_no_use_specs == 1)
  491.         int_error("Need 2 or 3 columns for polar data", this_plot->token);
  492.     }
  493.  
  494.     this_plot->num_iso_read = 0;
  495.     this_plot->has_grid_topology = TRUE;
  496.  
  497.     /* we ought to keep old memory - most likely case
  498.      * is a replot, so it will probably exactly fit into
  499.      * memory already allocated ?
  500.      */
  501.     if (this_plot->iso_crvs != NULL) {
  502.     struct iso_curve *icrv, *icrvs = this_plot->iso_crvs;
  503.  
  504.     while (icrvs) {
  505.         icrv = icrvs;
  506.         icrvs = icrvs->next;
  507.         iso_free(icrv);
  508.     }
  509.     this_plot->iso_crvs = NULL;
  510.     }
  511.     /* data file is already open */
  512.  
  513.     if (df_matrix)
  514.     xdatum = df_3dmatrix(this_plot);
  515.     else {
  516.     /*{{{  read surface from text file */
  517.     struct iso_curve *this_iso = iso_alloc(samples);
  518.     struct coordinate GPHUGE *cp;
  519.     double x, y, z;
  520.  
  521.     while ((j = df_readline(v, 3)) != DF_EOF) {
  522.         if (j == DF_SECOND_BLANK)
  523.         break;        /* two blank lines */
  524.         if (j == DF_FIRST_BLANK) {
  525.         /* one blank line */
  526.         if (pt_in_iso_crv == 0) {
  527.             if (xdatum == 0)
  528.             continue;
  529.             pt_in_iso_crv = xdatum;
  530.         }
  531.         if (xdatum > 0) {
  532.             this_iso->p_count = xdatum;
  533.             this_iso->next = this_plot->iso_crvs;
  534.             this_plot->iso_crvs = this_iso;
  535.             this_plot->num_iso_read++;
  536.  
  537.             if (xdatum != pt_in_iso_crv)
  538.             this_plot->has_grid_topology = FALSE;
  539.  
  540.             this_iso = iso_alloc(pt_in_iso_crv);
  541.             xdatum = 0;
  542.             ydatum++;
  543.         }
  544.         continue;
  545.         }
  546.         /* its a data point or undefined */
  547.  
  548.         if (xdatum >= this_iso->p_max) {
  549.         /*
  550.          * overflow about to occur. Extend size of points[] array. We
  551.          * either double the size, or add 1000 points, whichever is a
  552.          * smaller increment. Note i = p_max.
  553.          */
  554.         iso_extend(this_iso,
  555.                xdatum + (xdatum < 1000 ? xdatum : 1000));
  556.         }
  557.         cp = this_iso->points + xdatum;
  558.  
  559.         if (j == DF_UNDEFINED) {
  560.         cp->type = UNDEFINED;
  561.         continue;
  562.         }
  563.         cp->type = INRANGE;    /* unless we find out different */
  564.  
  565.         switch (mapping3d) {
  566.         case MAP3D_CARTESIAN:
  567.         switch (j) {
  568.         case 1:
  569.             x = xdatum;
  570.             y = ydatum;
  571.             z = v[0];
  572.             break;
  573.         case 3:
  574.             x = v[0];
  575.             y = v[1];
  576.             z = v[2];
  577.             break;
  578.         default:
  579.             {
  580.             char msg[80];
  581.             sprintf(msg, "Need 1 or 3 columns - line %d", df_line_number);
  582.             int_error(msg, this_plot->token);
  583.             return;    /* avoid gcc -Wuninitialised for x,y,z */
  584.             }
  585.         }
  586.         break;
  587.         case MAP3D_SPHERICAL:
  588.         if (j < 2)
  589.             int_error("Need 2 or 3 columns", this_plot->token);
  590.         if (j < 3)
  591.             v[2] = 1;    /* default radius */
  592.         if (angles_format == ANGLES_DEGREES) {
  593.             v[0] *= DEG2RAD;    /* Convert to radians. */
  594.             v[1] *= DEG2RAD;
  595.         }
  596.         x = v[2] * cos(v[0]) * cos(v[1]);
  597.         y = v[2] * sin(v[0]) * cos(v[1]);
  598.         z = v[2] * sin(v[1]);
  599.         break;
  600.         case MAP3D_CYLINDRICAL:
  601.         if (j < 2)
  602.             int_error("Need 2 or 3 columns", this_plot->token);
  603.         if (j < 3)
  604.             v[2] = 1;    /* default radius */
  605.         if (angles_format == ANGLES_DEGREES) {
  606.             v[0] *= DEG2RAD;    /* Convert to radians. */
  607.         }
  608.         x = v[2] * cos(v[0]);
  609.         y = v[2] * sin(v[0]);
  610.         z = v[1];
  611.         break;
  612.         default:
  613.         int_error("Internal error : Unknown mapping type", NO_CARET);
  614.         return;
  615.         }
  616.  
  617.         /* adjust for logscales. Set min/max and point types.
  618.          * store in cp
  619.          */
  620.         cp->type = INRANGE;
  621.         /* cannot use continue, as macro is wrapped in a loop.
  622.          * I regard this as correct goto use
  623.          */
  624.         STORE_WITH_LOG_AND_FIXUP_RANGE(cp->x, x, cp->type, x_axis, NOOP, goto come_here_if_undefined);
  625.         STORE_WITH_LOG_AND_FIXUP_RANGE(cp->y, y, cp->type, y_axis, NOOP, goto come_here_if_undefined);
  626.         STORE_WITH_LOG_AND_FIXUP_RANGE(cp->z, z, cp->type, z_axis, NOOP, goto come_here_if_undefined);
  627.  
  628.         /* some may complain, but I regard this as the correct use
  629.          * of goto
  630.          */
  631.       come_here_if_undefined:
  632.         ++xdatum;
  633.     }            /* end of whileloop - end of surface */
  634.  
  635.     if (xdatum > 0) {
  636.         this_plot->num_iso_read++;    /* Update last iso. */
  637.         this_iso->p_count = xdatum;
  638.  
  639.         this_iso->next = this_plot->iso_crvs;
  640.         this_plot->iso_crvs = this_iso;
  641.  
  642.         if (xdatum != pt_in_iso_crv)
  643.         this_plot->has_grid_topology = FALSE;
  644.  
  645.     } else {
  646.         iso_free(this_iso);    /* Free last allocation. */
  647.     }
  648.  
  649.     if (dgrid3d && this_plot->num_iso_read > 0)
  650.         grid_nongrid_data(this_plot);
  651.     /*}}} */
  652.     }
  653.  
  654.     if (this_plot->num_iso_read <= 1)
  655.     this_plot->has_grid_topology = FALSE;
  656.     if (this_plot->has_grid_topology && !hidden3d) {
  657.     struct iso_curve *new_icrvs = NULL;
  658.     int num_new_iso = this_plot->iso_crvs->p_count, len_new_iso = this_plot->num_iso_read;
  659.  
  660.     /* Now we need to set the other direction (pseudo) isolines. */
  661.     for (i = 0; i < num_new_iso; i++) {
  662.         struct iso_curve *new_icrv = iso_alloc(len_new_iso);
  663.  
  664.         new_icrv->p_count = len_new_iso;
  665.  
  666.         for (j = 0, this_iso = this_plot->iso_crvs;
  667.          this_iso != NULL;
  668.          j++, this_iso = this_iso->next) {
  669.         /* copy whole point struct to get type too.
  670.          * wasteful for windows, with padding */
  671.         /* more efficient would be extra pointer to same struct */
  672.         new_icrv->points[j] = this_iso->points[i];
  673.         }
  674.  
  675.         new_icrv->next = new_icrvs;
  676.         new_icrvs = new_icrv;
  677.     }
  678.  
  679.     /* Append the new iso curves after the read ones. */
  680.     for (this_iso = this_plot->iso_crvs;
  681.          this_iso->next != NULL;
  682.          this_iso = this_iso->next);
  683.     this_iso->next = new_icrvs;
  684.     }
  685. }
  686.  
  687.  
  688.  
  689. static void print_3dtable(pcount)
  690. int pcount;
  691. {
  692.     register struct surface_points *this_plot;
  693.     int i, curve, surface;
  694.     struct iso_curve *icrvs;
  695.     struct coordinate GPHUGE *points;
  696.  
  697.     for (surface = 0, this_plot = first_3dplot; surface < pcount;
  698.      this_plot = this_plot->next_sp, surface++) {
  699.     fprintf(gpoutfile, "\n#Surface %d of %d surfaces\n", surface, pcount);
  700.     icrvs = this_plot->iso_crvs;
  701.     curve = 0;
  702.  
  703.     if (draw_surface) {
  704.         /* only the curves in one direction */
  705.         while (icrvs && curve < this_plot->num_iso_read) {
  706.         fprintf(gpoutfile, "\n#IsoCurve %d, %d points\n#x y z type\n",
  707.             curve, icrvs->p_count);
  708.         for (i = 0, points = icrvs->points; i < icrvs->p_count; i++) {
  709.             fprintf(gpoutfile, "%g %g %g %c\n",
  710.                 points[i].x,
  711.                 points[i].y,
  712.                 points[i].z,
  713.                 points[i].type == INRANGE ? 'i'
  714.                 : points[i].type == OUTRANGE ? 'o'
  715.                 : 'u');
  716.         }
  717.         icrvs = icrvs->next;
  718.         curve++;
  719.         }
  720.         putc('\n', gpoutfile);
  721.     }
  722.     if (draw_contour) {
  723.         int number = 0;
  724.         struct gnuplot_contours *c = this_plot->contours;
  725.         while (c) {
  726.         int count = c->num_pts;
  727.         struct coordinate GPHUGE *p = c->coords;
  728.         if (c->isNewLevel)
  729.             /* dont display count - contour split across chunks */
  730.             /* put # in case user wants to use it for a plot */
  731.             /* double blank line to allow plot ... index ... */
  732.             fprintf(gpoutfile, "\n# Contour %d, label: %s\n", number++, c->label);
  733.         for (; --count >= 0; ++p)
  734.             fprintf(gpoutfile, "%g %g %g\n", p->x, p->y, p->z);
  735.         /* blank line between segments of same contour */
  736.         putc('\n', gpoutfile);
  737.         c = c->next;
  738.         }
  739.     }
  740.     }
  741.     fflush(gpoutfile);
  742. }
  743.  
  744.  
  745.  
  746. #define SET_DUMMY_RANGE(AXIS) \
  747. do{\
  748.  if (parametric || polar) { \
  749.   t_min = tmin; t_max = tmax;\
  750.  } else if (log_array[AXIS]) {\
  751.   if (min_array[AXIS] <= 0.0 || max_array[AXIS] <= 0.0)\
  752.    int_error("x/x2 range must be greater than 0 for log scale!", NO_CARET);\
  753.   t_min = log(min_array[AXIS])/log_base_array[AXIS]; t_max = log(max_array[AXIS])/log_base_array[AXIS];\
  754.  } else {\
  755.   t_min = min_array[AXIS]; t_max = max_array[AXIS];\
  756.  }\
  757.  t_step = (t_max - t_min) / (samples - 1); \
  758. }while(0)
  759.  
  760.  
  761. /*
  762.  * This parses the splot command after any range specifications. To support
  763.  * autoscaling on the x/z axis, we want any data files to define the x/y
  764.  * range, then to plot any functions using that range. We thus parse the
  765.  * input twice, once to pick up the data files, and again to pick up the
  766.  * functions. Definitions are processed twice, but that won't hurt.
  767.  * div - okay, it doesn't hurt, but every time an option as added for
  768.  * datafiles, code to parse it has to be added here. Change so that
  769.  * we store starting-token in the plot structure.
  770.  */
  771. static void eval_3dplots()
  772. {
  773.     int i, j;
  774.     struct surface_points **tp_3d_ptr;
  775.     int start_token, end_token;
  776.     int begin_token;
  777.     TBOOLEAN some_data_files = FALSE, some_functions = FALSE;
  778.     int plot_num, line_num, point_num, crnt_param = 0;    /* 0 = z, 1 = y, 2 = x */
  779.     char *xtitle;
  780.     char *ytitle;
  781.  
  782.     /* Reset first_3dplot. This is usually done at the end of this function.
  783.      * If there is an error within this function, the memory is left allocated,
  784.      * since we cannot call sp_free if the list is incomplete
  785.      */
  786.     first_3dplot = NULL;
  787.  
  788.     /* put stuff into arrays to simplify access */
  789.     INIT_ARRAYS(FIRST_X_AXIS, xmin, xmax, autoscale_lx, is_log_x, base_log_x, log_base_log_x, 0);
  790.     INIT_ARRAYS(FIRST_Y_AXIS, ymin, ymax, autoscale_ly, is_log_y, base_log_y, log_base_log_y, 0);
  791.     INIT_ARRAYS(FIRST_Z_AXIS, zmin, zmax, autoscale_lz, is_log_z, base_log_z, log_base_log_z, 1);
  792.  
  793.     x_axis = FIRST_X_AXIS;
  794.     y_axis = FIRST_Y_AXIS;
  795.     z_axis = FIRST_Z_AXIS;
  796.  
  797.     tp_3d_ptr = &(first_3dplot);
  798.     plot_num = 0;
  799.     line_num = 0;        /* default line type */
  800.     point_num = 0;        /* default point type */
  801.  
  802.     xtitle = NULL;
  803.     ytitle = NULL;
  804.  
  805.     begin_token = c_token;
  806.  
  807. /*** First Pass: Read through data files ***/
  808.     /*
  809.      * This pass serves to set the x/yranges and to parse the command, as
  810.      * well as filling in every thing except the function data. That is done
  811.      * after the x/yrange is defined.
  812.      */
  813.     while (TRUE) {
  814.     if (END_OF_COMMAND)
  815.         int_error("function to plt3d expected", c_token);
  816.  
  817.     start_token = c_token;
  818.  
  819.     if (is_definition(c_token)) {
  820.         define();
  821.     } else {
  822.         int specs;
  823.         struct surface_points *this_plot;
  824.  
  825.         if (isstring(c_token)) {    /* data file to plot */
  826.  
  827.         /*{{{  data file */
  828.         if (parametric && crnt_param != 0)
  829.             int_error("previous parametric function not fully specified", c_token);
  830.  
  831.         if (!some_data_files) {
  832.             if (autoscale_lx & 1) {
  833.             min_array[FIRST_X_AXIS] = VERYLARGE;
  834.             }
  835.             if (autoscale_lx & 2) {
  836.             max_array[FIRST_X_AXIS] = -VERYLARGE;
  837.             }
  838.             if (autoscale_ly & 1) {
  839.             min_array[FIRST_Y_AXIS] = VERYLARGE;
  840.             }
  841.             if (autoscale_ly & 2) {
  842.             max_array[FIRST_Y_AXIS] = -VERYLARGE;
  843.             }
  844.             some_data_files = TRUE;
  845.         }
  846.         if (*tp_3d_ptr)
  847.             this_plot = *tp_3d_ptr;
  848.         else {        /* no memory malloc()'d there yet */
  849.             /* Allocate enough isosamples and samples */
  850.             this_plot = sp_alloc(0, 0, 0, 0);
  851.             *tp_3d_ptr = this_plot;
  852.         }
  853.  
  854.         this_plot->plot_type = DATA3D;
  855.         this_plot->plot_style = data_style;
  856.  
  857.         specs = df_open(3);
  858.         /* parses all datafile-specific modifiers */
  859.         /* we will load the data after parsing title,with,... */
  860.  
  861.         /* for capture to key */
  862.         this_plot->token = end_token = c_token - 1;
  863.  
  864.         /* this_plot->token is temporary, for errors in get_3ddata() */
  865.  
  866.         if (datatype[FIRST_X_AXIS] == TIME) {
  867.             if (specs < 3)
  868.             int_error("Need full using spec for x time data",
  869.                   c_token);
  870.             df_timecol[0] = 1;
  871.         }
  872.         if (datatype[FIRST_Y_AXIS] == TIME) {
  873.             if (specs < 3)
  874.             int_error("Need full using spec for y time data",
  875.                   c_token);
  876.             df_timecol[1] = 1;
  877.         }
  878.         if (datatype[FIRST_Z_AXIS] == TIME) {
  879.             if (specs < 3)
  880.             df_timecol[0] = 1;
  881.             else
  882.             df_timecol[2] = 1;
  883.         }
  884.         /*}}} */
  885.  
  886.         } else {        /* function to plot */
  887.  
  888.         /*{{{  function */
  889.         ++plot_num;
  890.         if (parametric)    {
  891.             /* Rotate between x/y/z axes */
  892.             /* +2 same as -1, but beats -ve problem */
  893.             crnt_param = (crnt_param + 2) % 3;
  894.         }
  895.  
  896.         if (*tp_3d_ptr) {
  897.             this_plot = *tp_3d_ptr;
  898.             if (!hidden3d)
  899.             sp_replace(this_plot, samples_1, iso_samples_1,
  900.                    samples_2, iso_samples_2);
  901.             else
  902.             sp_replace(this_plot, iso_samples_1, 0,
  903.                    0, iso_samples_2);
  904.  
  905.         } else {    /* no memory malloc()'d there yet */
  906.             /* Allocate enough isosamples and samples */
  907.             if (!hidden3d)
  908.             this_plot = sp_alloc(samples_1, iso_samples_1,
  909.                          samples_2, iso_samples_2);
  910.             else
  911.             this_plot = sp_alloc(iso_samples_1, 0,
  912.                          0, iso_samples_2);
  913.             *tp_3d_ptr = this_plot;
  914.         }
  915.  
  916.         this_plot->plot_type = FUNC3D;
  917.         this_plot->has_grid_topology = TRUE;
  918.         this_plot->plot_style = func_style;
  919.         this_plot->num_iso_read = iso_samples_2;
  920.         dummy_func = &plot_func;
  921.         plot_func.at = temp_at();
  922.         dummy_func = NULL;
  923.         /* ignore it for now */
  924.         some_functions = TRUE;
  925.         end_token = c_token - 1;
  926.         /*}}} */
  927.  
  928.         }            /* end of IS THIS A FILE OR A FUNC block */
  929.  
  930.  
  931.         /*{{{  title */
  932.         if (this_plot->title) {
  933.         free(this_plot->title);
  934.         this_plot->title = NULL;
  935.         }
  936.         if (almost_equals(c_token, "t$itle")) {
  937.         /*                      if (!isstring(++c_token))
  938.            int_error("Expected title", c_token);
  939.            m_quote_capture(&(this_plot->title), c_token, c_token);
  940.          */
  941.         if (parametric) {
  942.             if (crnt_param != 0)
  943.             int_error("\"title\" allowed only after parametric function fully specified", c_token);
  944.             else {
  945.             if (xtitle != NULL)
  946.                 xtitle[0] = NUL;    /* Remove default title . */
  947.             if (ytitle != NULL)
  948.                 ytitle[0] = NUL;    /* Remove default title . */
  949.             }
  950.         }
  951.         if (isstring(++c_token))
  952.             m_quote_capture(&(this_plot->title), c_token, c_token);
  953.         else
  954.             int_error("expecting \"title\" for plot", c_token);
  955.         /* end of new method */
  956.         ++c_token;
  957.         } else if (almost_equals(c_token, "not$itle")) {
  958.         if (xtitle != NULL)
  959.             xtitle[0] = '\0';
  960.         if (ytitle != NULL)
  961.             ytitle[0] = '\0';
  962.         /*   this_plot->title = NULL;   */
  963.         ++c_token;
  964.         } else {
  965.         m_capture(&(this_plot->title), start_token, end_token);
  966.         if (crnt_param == 2)
  967.             xtitle = this_plot->title;
  968.         else if (crnt_param == 1)
  969.             ytitle = this_plot->title;
  970.         }
  971.         /*}}} */
  972.  
  973.  
  974.         /*{{{  line types, widths, ... */
  975.         this_plot->lp_properties.l_type = line_num;
  976.         this_plot->lp_properties.p_type = point_num;
  977.  
  978.         if (almost_equals(c_token, "w$ith")) {
  979.         this_plot->plot_style = get_style();
  980.         }
  981.         /* pick up line/point specs
  982.          * - point spec allowed if style uses points, ie style&2 != 0
  983.          * - keywords are optional
  984.          */
  985.         LP_PARSE(this_plot->lp_properties, 1, this_plot->plot_style & 2,
  986.              line_num, point_num);
  987.  
  988.         /* allow old-style syntax too - ignore case lt 3 4 for example */
  989.         if (!equals(c_token, ",") && !END_OF_COMMAND) {
  990.         struct value t;
  991.         this_plot->lp_properties.l_type =
  992.             this_plot->lp_properties.p_type = (int) real(const_express(&t)) - 1;
  993.  
  994.         if (!equals(c_token, ",") && !END_OF_COMMAND)
  995.             this_plot->lp_properties.p_type = (int) real(const_express(&t)) - 1;
  996.         }
  997.         if (this_plot->plot_style & 2)    /* lines, linesp, ... */
  998.         if (crnt_param == 0)
  999.             point_num +=
  1000.             1 + (draw_contour != 0)
  1001.             + (hidden3d != 0);
  1002.  
  1003.         if (crnt_param == 0)
  1004.         line_num += 1 + (draw_contour != 0)
  1005.             + (hidden3d != 0);
  1006.         /*}}} */
  1007.  
  1008.  
  1009.         /* now get the data... having to think hard here...
  1010.          * first time through, we fill in this_plot. For second
  1011.          * surface in file, we have to allocate another surface
  1012.          * struct. BUT we may allocate this store only to
  1013.          * find that it is merely some blank lines at end of file
  1014.          * tp_3d_ptr is still pointing at next field of prev. plot,
  1015.          * before :    prev_or_first -> this_plot -> possible_preallocated_store
  1016.          *                tp_3d_ptr--^
  1017.          * after  :    prev_or_first -> first -> second -> last -> possibly_more_store
  1018.          *                                        tp_3d_ptr ----^
  1019.          * if file is empty, tp_3d_ptr is not moved. this_plot continues
  1020.          * to point at allocated storage, but that will be reused later
  1021.          */
  1022.  
  1023.         assert(this_plot == *tp_3d_ptr);
  1024.  
  1025.         if (this_plot->plot_type == DATA3D) {
  1026.         /*{{{  read data */
  1027.         /* remember settings for second surface in file */
  1028.         struct lp_style_type *these_props = &(this_plot->lp_properties);
  1029.         enum PLOT_STYLE this_style = this_plot->plot_style;
  1030.  
  1031.         int this_token = this_plot->token;
  1032.         while (!df_eof) {
  1033.             this_plot = *tp_3d_ptr;
  1034.             assert(this_plot != NULL);
  1035.  
  1036.             /* dont move tp_3d_ptr until we are sure we
  1037.              * have read a surface
  1038.              */
  1039.             
  1040.             /* used by get_3ddata() */
  1041.             this_plot->token = this_token;
  1042.  
  1043.             get_3ddata(this_plot);
  1044.             /* for second pass */
  1045.             this_plot->token = c_token;
  1046.  
  1047.             if (this_plot->num_iso_read == 0)
  1048.             /* probably df_eof, in which case we
  1049.              * will leave loop. if not eof, then
  1050.              * how come we got no surface ? - retry
  1051.              * in neither case do we update tp_3d_ptr
  1052.              */
  1053.             continue;
  1054.  
  1055.             /* okay, we have read a surface */
  1056.             ++plot_num;
  1057.             tp_3d_ptr = &(this_plot->next_sp);
  1058.             if (df_eof)
  1059.             break;
  1060.  
  1061.             /* there might be another surface so allocate
  1062.              * and prepare another surface structure
  1063.              * This does no harm if in fact there are
  1064.              * no more surfaces to read
  1065.              */
  1066.  
  1067.             if ((this_plot = *tp_3d_ptr) != NULL) {
  1068.             if (this_plot->title) {
  1069.                 free(this_plot->title);
  1070.                 this_plot->title = NULL;
  1071.             }
  1072.             } else {
  1073.             /* Allocate enough isosamples and samples */
  1074.             this_plot = *tp_3d_ptr = sp_alloc(0, 0, 0, 0);
  1075.             }
  1076.  
  1077.             this_plot->plot_type = DATA3D;
  1078.             this_plot->plot_style = this_style;
  1079.             /* Struct copy */
  1080.             this_plot->lp_properties = *these_props;
  1081.         }
  1082.         df_close();
  1083.         /*}}} */
  1084.  
  1085.         } else {        /* not a data file */
  1086.         tp_3d_ptr = &(this_plot->next_sp);
  1087.         this_plot->øÄ
  1088.     /* store for second pass */
  1089.         }
  1090.  
  1091.     }        /* !is_definition() : end of scope of this_plot */
  1092.  
  1093.  
  1094.     if (equals(c_token, ","))
  1095.         c_token++;
  1096.     else
  1097.         break;
  1098.  
  1099.     }                /* while(TRUE), ie first pass */
  1100.  
  1101.  
  1102.     if (parametric && crnt_param != 0)
  1103.     int_error("parametric function not fully specified", NO_CARET);
  1104.  
  1105.  
  1106. /*** Second Pass: Evaluate the functions ***/
  1107.     /*
  1108.      * Everything is defined now, except the function data. We expect no
  1109.      * syntax errors, etc, since the above parsed it all. This makes the code
  1110.      * below simpler. If autoscale_ly, the yrange may still change.
  1111.      * - eh ?  - z can still change.  x/y/z can change if we are parametric ??
  1112.      */
  1113.  
  1114.     if (some_functions) {
  1115.  
  1116.     /* I've changed the controlled variable in fn plots to u_min etc since
  1117.      * it's easier for me to think parametric - 'normal' plot is after all
  1118.      * a special case. I was confused about x_min being both minimum of
  1119.      * x values found, and starting value for fn plots.
  1120.      */
  1121.     register double u_min, u_max, u_step, v_min, v_max, v_step;
  1122.     double uisodiff, visodiff;
  1123.     struct surface_points *this_plot;
  1124.  
  1125.  
  1126.     if (!parametric) {
  1127.         /*{{{  check ranges */
  1128.         /* give error if xrange badly set from missing datafile error
  1129.          * parametric fn can still set ranges
  1130.          * if there are no fns, we'll report it later as 'nothing to plot'
  1131.          */
  1132.  
  1133.         if (min_array[FIRST_X_AXIS] == VERYLARGE || 
  1134.         max_array[FIRST_X_AXIS] == -VERYLARGE) {
  1135.         int_error("x range is invalid", c_token);
  1136.         }
  1137.         if (min_array[FIRST_Y_AXIS] == VERYLARGE || 
  1138.         max_array[FIRST_Y_AXIS] == -VERYLARGE) {
  1139.         int_error("y range is invalid", c_token);
  1140.         }
  1141.         /* check that xmin -> xmax is not too small */
  1142.         fixup_range(FIRST_X_AXIS, "x");
  1143.         fixup_range(FIRST_Y_AXIS, "y");
  1144.         /*}}} */
  1145.     }
  1146.     if (parametric && !some_data_files) {
  1147.         /*{{{  set ranges */
  1148.         /* parametric fn can still change x/y range */
  1149.         if (autoscale_lx & 1)
  1150.         min_array[FIRST_X_AXIS] = VERYLARGE;
  1151.         if (autoscale_lx & 2)
  1152.         max_array[FIRST_X_AXIS] = -VERYLARGE;
  1153.         if (autoscale_ly & 1)
  1154.         min_array[FIRST_Y_AXIS] = VERYLARGE;
  1155.         if (autoscale_ly & 2)
  1156.         max_array[FIRST_Y_AXIS] = -VERYLARGE;
  1157.         /*}}} */
  1158.     }
  1159.     if (parametric) {
  1160.         u_min = umin;
  1161.         u_max = umax;
  1162.         v_min = vmin;
  1163.         v_max = vmax;
  1164.     } else {
  1165.         /*{{{  figure ranges, taking logs etc into account */
  1166.         if (is_log_x) {
  1167.         if (min_array[FIRST_X_AXIS] <= 0.0 ||
  1168.             max_array[FIRST_X_AXIS] <= 0.0)
  1169.             int_error("x range must be greater than 0 for log scale!",
  1170.                   NO_CARET);
  1171.         u_min = log(min_array[FIRST_X_AXIS]) / log_base_log_x;
  1172.         u_max = log(max_array[FIRST_X_AXIS]) / log_base_log_x;
  1173.         } else {
  1174.         u_min = min_array[FIRST_X_AXIS];
  1175.         u_max = max_array[FIRST_X_AXIS];
  1176.         }
  1177.  
  1178.         if (is_log_y) {
  1179.         if (min_array[FIRST_Y_AXIS] <= 0.0 ||
  1180.             max_array[FIRST_Y_AXIS] <= 0.0) {
  1181.             int_error("y range must be greater than 0 for log scale!",
  1182.                   NO_CARET);
  1183.         }
  1184.         v_min = log(min_array[FIRST_Y_AXIS]) / log_base_log_y;
  1185.         v_max = log(max_array[FIRST_Y_AXIS]) / log_base_log_y;
  1186.         } else {
  1187.         v_min = min_array[FIRST_Y_AXIS];
  1188.         v_max = max_array[FIRST_Y_AXIS];
  1189.         }
  1190.         /*}}} */
  1191.     }
  1192.  
  1193.  
  1194.     if (samples_1 < 2 || samples_2 < 2 || iso_samples_1 < 2 ||
  1195.         iso_samples_2 < 2) {
  1196.         int_error("samples or iso_samples < 2. Must be at least 2.",
  1197.               NO_CARET);
  1198.     }
  1199.  
  1200.     /* start over */
  1201.     this_plot = first_3dplot;
  1202.     c_token = begin_token;
  1203.  
  1204.     /* why do attributes of this_plot matter ? */
  1205.  
  1206.     if (this_plot && this_plot->has_grid_topology && hidden3d) {
  1207.         u_step = (u_max - u_min) / (iso_samples_1 - 1);
  1208.         v_step = (v_max - v_min) / (iso_samples_2 - 1);
  1209.     } else {
  1210.         u_step = (u_max - u_min) / (samples_1 - 1);
  1211.         v_step = (v_max - v_min) / (samples_2 - 1);
  1212.     }
  1213.  
  1214.     uisodiff = (u_max - u_min) / (iso_samples_1 - 1);
  1215.     visodiff = (v_max - v_min) / (iso_samples_2 - 1);
  1216.  
  1217.  
  1218.     /* Read through functions */
  1219.     while (TRUE) {
  1220.         if (is_definition(c_token)) {
  1221.         define();
  1222.         } else {
  1223.  
  1224.         if (!isstring(c_token)) {    /* func to plot */
  1225.             /*{{{  evaluate function */
  1226.             struct iso_curve *this_iso = this_plot->iso_crvs;
  1227.             struct coordinate GPHUGE *points = this_iso->points;
  1228.             int num_sam_to_use, num_iso_to_use;
  1229.  
  1230.             if (parametric)
  1231.             crnt_param = (crnt_param + 2) % 3;
  1232.  
  1233.             dummy_func = &plot_func;
  1234.             plot_func.at = temp_at();    /* reparse function */
  1235.             dummy_func = NULL;
  1236.             num_iso_to_use = iso_samples_2;
  1237.  
  1238.             if (!(this_plot->has_grid_topology && hidden3d))
  1239.             num_sam_to_use = samples_1;
  1240.             else
  1241.             num_sam_to_use = iso_samples_1;
  1242.  
  1243.             for (j = 0; j < num_iso_to_use; j++) {
  1244.             double y = v_min + j * visodiff;
  1245.             /* if (is_log_y) PEM fix logscale y axis */
  1246.             /* y = pow(log_base_log_y,y); 26-Sep-89 */
  1247.             /* parametric => NOT a log quantity (?) */
  1248.             (void) Gcomplex(&plot_func.dummy_values[1],
  1249.             !parametric && is_log_y ? pow(base_log_y, y) : y,
  1250.                     0.0);
  1251.  
  1252.             for (i = 0; i < num_sam_to_use; i++) {
  1253.                 double x = u_min + i * u_step;
  1254.                 struct value a;
  1255.                 double temp;
  1256.  
  1257.                 /* if (is_log_x) PEM fix logscale x axis */
  1258.                 /* x = pow(base_log_x,x); 26-Sep-89 */
  1259.                 /* parametric => NOT a log quantity (?) */
  1260.                 (void) Gcomplex(&plot_func.dummy_values[0],
  1261.                         !parametric && is_log_x ? pow(base_log_x, x) : x, 0.0);
  1262.  
  1263.                 points[i].x = x;
  1264.                 points[i].y = y;
  1265.  
  1266.                 evaluate_at(plot_func.at, &a);
  1267.  
  1268.                 if (undefined || (fabs(imag(&a)) > zero)) {
  1269.                 points[i].type = UNDEFINED;
  1270.                 continue;
  1271.                 }
  1272.                 temp = real(&a);
  1273.  
  1274.                 points[i].type = INRANGE;
  1275.                 STORE_WITH_LOG_AND_FIXUP_RANGE(points[i].z, temp, points[i].type, crnt_param, NOOP, NOOP);
  1276.  
  1277.             }
  1278.             this_iso->p_count = num_sam_to_use;
  1279.             this_iso = this_iso->next;
  1280.             points = this_iso ? this_iso->points : NULL;
  1281.             }
  1282.  
  1283.             if (!(this_plot->has_grid_topology && hidden3d)) {
  1284.             num_iso_to_use = iso_samples_1;
  1285.             num_sam_to_use = samples_2;
  1286.             for (i = 0; i < num_iso_to_use; i++) {
  1287.                 double x = u_min + i * uisodiff;
  1288.                 /* if (is_log_x) PEM fix logscale x axis */
  1289.                 /* x = pow(base_log_x,x); 26-Sep-89 */
  1290.                 /* if parametric, no logs involved - 3.6 */
  1291.                 (void) Gcomplex(&plot_func.dummy_values[0],
  1292.                         (!parametric && is_log_x) ? pow(base_log_x, x) : x, 0.0);
  1293.  
  1294.                 for (j = 0; j < num_sam_to_use; j++) {
  1295.                 double y = v_min + j * v_step;
  1296.                 struct value a;
  1297.                 double temp;
  1298.                 /* if (is_log_y) PEM fix logscale y axis */
  1299.                 /* y = pow(base_log_y,y); 26-Sep-89 */
  1300.                 (void) Gcomplex(&plot_func.dummy_values[1],
  1301.                         (!parametric && is_log_y) ? pow(base_log_y, y) : y, 0.0);
  1302.  
  1303.                 points[j].x = x;
  1304.                 points[j].y = y;
  1305.  
  1306.                 evaluate_at(plot_func.at, &a);
  1307.  
  1308.                 if (undefined || (fabs(imag(&a)) > zero)) {
  1309.                     points[j].type = UNDEFINED;
  1310.                     continue;
  1311.                 }
  1312.                 temp = real(&a);
  1313.  
  1314.                 points[j].type = INRANGE;
  1315.                 STORE_WITH_LOG_AND_FIXUP_RANGE(points[j].z, temp, points[j].type, crnt_param, NOOP, NOOP);
  1316.                 }
  1317.                 this_iso->p_count = num_sam_to_use;
  1318.                 this_iso = this_iso->next;
  1319.                 points = this_iso ? this_iso->points : NULL;
  1320.             }
  1321.             }
  1322.             /*}}} */
  1323.         }        /* end of ITS A FUNCTION TO PLOT */
  1324.  
  1325.         /* we saved it from first pass */ 
  1326.         c_token = this_plot->token;
  1327.  
  1328.         /* one data file can make several plots */
  1329.         do
  1330.             this_plot = this_plot->next_sp;
  1331.         while (this_plot && this_plot->token == c_token);
  1332.  
  1333.         }            /* !is_definition */
  1334.  
  1335.         if (equals(c_token, ","))
  1336.         c_token++;
  1337.         else
  1338.         break;
  1339.  
  1340.     }            /* while(TRUE) */
  1341.  
  1342.  
  1343.     if (parametric) {
  1344.         /* Now actually fix the plot triplets to be single plots. */
  1345.         parametric_3dfixup(first_3dplot, &plot_num);
  1346.     }
  1347.     }                /* some functions */
  1348.     /* if first_3dplot is NULL, we have no functions or data at all.
  1349.        * This can happen, if you type "splot x=5", since x=5 is a
  1350.        * variable assignment
  1351.      */
  1352.     if (plot_num == 0 || first_3dplot == NULL) {
  1353.     int_error("no functions or data to plot", c_token);
  1354.     }
  1355.     if (min_array[FIRST_X_AXIS] == VERYLARGE ||
  1356.     max_array[FIRST_X_AXIS] == -VERYLARGE ||
  1357.     min_array[FIRST_Y_AXIS] == VERYLARGE ||
  1358.     max_array[FIRST_Y_AXIS] == -VERYLARGE ||
  1359.     min_array[FIRST_Z_AXIS] == VERYLARGE ||
  1360.     max_array[FIRST_Z_AXIS] == -VERYLARGE)
  1361.     int_error("All points undefined", NO_CARET);
  1362.  
  1363.     fixup_range(FIRST_X_AXIS, "x");
  1364.     fixup_range(FIRST_Y_AXIS, "y");
  1365.     fixup_range(FIRST_Z_AXIS, "z");
  1366.  
  1367.     FIXUP_RANGE_FOR_LOG(FIRST_X_AXIS, x);
  1368.     FIXUP_RANGE_FOR_LOG(FIRST_Y_AXIS, y);
  1369.     FIXUP_RANGE_FOR_LOG(FIRST_Z_AXIS, z);
  1370.  
  1371.     /* last parameter should take plot size into effect...
  1372.      * probably needs to be moved to graph3d.c
  1373.      * in the meantime, a value of 20 gives same behaviour
  1374.      * as 3.5 which will do for the moment
  1375.      */
  1376.  
  1377.     if (xtics)
  1378.     setup_tics(FIRST_X_AXIS, &xticdef, xformat, 20);
  1379.     if (ytics)
  1380.     setup_tics(FIRST_Y_AXIS, &yticdef, yformat, 20);
  1381.     if (ztics)
  1382.     setup_tics(FIRST_Z_AXIS, &zticdef, zformat, 20);
  1383.  
  1384. #define WRITEBACK(axis,min,max) \
  1385. if(range_flags[axis]&RANGE_WRITEBACK) \
  1386.   {if (auto_array[axis]&1) min = min_array[axis]; \
  1387.    if (auto_array[axis]&2) max = max_array[axis]; \
  1388.   }
  1389.  
  1390.     WRITEBACK(FIRST_X_AXIS, xmin, xmax);
  1391.     WRITEBACK(FIRST_Y_AXIS, ymin, ymax);
  1392.     WRITEBACK(FIRST_Z_AXIS, zmin, zmax);
  1393.  
  1394.     if (plot_num == 0 || first_3dplot == NULL) {
  1395.         int_error("no functions or data to plot", c_token);
  1396.     }
  1397.     /* Creates contours if contours are to be plotted as well. */
  1398.  
  1399.     if (draw_contour) {
  1400.     struct surface_points *this_plot;
  1401.     for (this_plot = first_3dplot, i = 0;
  1402.          i < plot_num;
  1403.          this_plot = this_plot->next_sp, i++) {
  1404.  
  1405.         if (this_plot->contours) {
  1406.         struct gnuplot_contours *cntrs = this_plot->contours;
  1407.  
  1408.         while (cntrs) {
  1409.             struct gnuplot_contours *cntr = cntrs;
  1410.             cntrs = cntrs->next;
  1411.             free(cntr->coords);
  1412.             free(cntr);
  1413.         }
  1414.         }
  1415.         /* Make sure this one can be contoured. */
  1416.         if (!this_plot->has_grid_topology) {
  1417.         this_plot->contours = NULL;
  1418.         fputs("Notice: cannot contour non grid data!\n", stderr);
  1419.         /* changed from int_error by recommendation of
  1420.          * rkc@xn.ll.mit.edu
  1421.          */
  1422.         } else if (this_plot->plot_type == DATA3D) {
  1423.         this_plot->contours = contour(
  1424.                          this_plot->num_iso_read,
  1425.                          this_plot->iso_crvs,
  1426.                          contour_levels, contour_pts,
  1427.                          contour_kind, contour_order,
  1428.                            levels_kind, levels_list);
  1429.         } else {
  1430.         this_plot->contours = contour(iso_samples_2,
  1431.                           this_plot->iso_crvs,
  1432.                           contour_levels, contour_pts,
  1433.                           contour_kind, contour_order,
  1434.                           levels_kind, levels_list);
  1435.         }
  1436.     }
  1437.     }                /* draw_contour */
  1438.     /* perform the plot */
  1439.     if (strcmp(term->name, "table") == 0)
  1440.     print_3dtable(plot_num);
  1441.     else {
  1442.     START_LEAK_CHECK();    /* assert no memory leaks here ! */
  1443.     do_3dplot(first_3dplot, plot_num);
  1444.     END_LEAK_CHECK();
  1445.     }
  1446.  
  1447.     /* if we get here, all went well, so record the line for replot */
  1448.     if (plot_token != -1) {
  1449.     /* note that m_capture also frees the old replot_line */
  1450.     m_capture(&replot_line, plot_token, c_token - 1);
  1451.     plot_token = -1;
  1452.     }
  1453.     sp_free(first_3dplot);
  1454.     first_3dplot = NULL;
  1455. }
  1456.  
  1457.  
  1458.  
  1459.  
  1460.  
  1461. static void parametric_3dfixup(start_plot, plot_num)
  1462. struct surface_points *start_plot;
  1463. int *plot_num;
  1464. /*
  1465.  * The hardest part of this routine is collapsing the FUNC plot types in the
  1466.  * list (which are gauranteed to occur in (x,y,z) triplets while preserving
  1467.  * the non-FUNC type plots intact.  This means we have to work our way
  1468.  * through various lists.  Examples (hand checked):
  1469.  * start_plot:F1->F2->F3->NULL ==> F3->NULL
  1470.  * start_plot:F1->F2->F3->F4->F5->F6->NULL ==> F3->F6->NULL
  1471.  * start_plot:F1->F2->F3->D1->D2->F4->F5->F6->D3->NULL ==>
  1472.  * F3->D1->D2->F6->D3->NULL
  1473.  */
  1474. {
  1475. /*
  1476.  * I initialized *free_list with NULL, because my compiler warns some lines
  1477.  * later that it might be uninited. The code however seems to not access that
  1478.  * line in that case, but if I'm right, my change is OK and if not, this is a
  1479.  * serious bug in the code.
  1480.  *
  1481.  * x and y ranges now fixed in eval_3dplots
  1482.  */
  1483.     struct surface_points *xp, *new_list, *free_list = NULL;
  1484.     struct surface_points **last_pointer = &new_list;
  1485.  
  1486.     int i, tlen, surface;
  1487.     char *new_title;
  1488.  
  1489.     /*
  1490.      * Ok, go through all the plots and move FUNC3D types together.  Note:
  1491.      * this originally was written to look for a NULL next pointer, but
  1492.      * gnuplot wants to be sticky in grabbing memory and the right number of
  1493.      * items in the plot list is controlled by the plot_num variable.
  1494.      * 
  1495.      * Since gnuplot wants to do this sticky business, a free_list of
  1496.      * surface_points is kept and then tagged onto the end of the plot list
  1497.      * as this seems more in the spirit of the original memory behavior than
  1498.      * simply freeing the memory.  I'm personally not convinced this sort of
  1499.      * concern is worth it since the time spent computing points seems to
  1500.      * dominate any garbage collecting that might be saved here...
  1501.      */
  1502.     new_list = xp = start_plot;
  1503.     for (surface = 0; surface < *plot_num; surface++) {
  1504.     if (xp->plot_type == FUNC3D) {
  1505.         struct surface_points *yp = xp->next_sp;
  1506.         struct surface_points *zp = yp->next_sp;
  1507.  
  1508.         /* Here's a FUNC3D parametric function defined as three parts.
  1509.          * Go through all the points and assign the x's and y's from xp and
  1510.          * yp to zp. min/max already done
  1511.          */
  1512.         struct iso_curve *xicrvs = xp->iso_crvs;
  1513.         struct iso_curve *yicrvs = yp->iso_crvs;
  1514.         struct iso_curve *zicrvs = zp->iso_crvs;
  1515.  
  1516.         (*plot_num) -= 2;
  1517.  
  1518.         assert(INRANGE < OUTRANGE && OUTRANGE < UNDEFINED);
  1519.  
  1520.         while (zicrvs) {
  1521.         struct coordinate GPHUGE *xpoints = xicrvs->points,
  1522.          GPHUGE * ypoints = yicrvs->points, GPHUGE * zpoints = zicrvs->points;
  1523.         for (i = 0; i < zicrvs->p_count; ++i) {
  1524.             zpoints[i].x = xpoints[i].z;
  1525.             zpoints[i].y = ypoints[i].z;
  1526.             if (zpoints[i].type < xpoints[i].type)
  1527.             zpoints[i].type = xpoints[i].type;
  1528.             if (zpoints[i].type < ypoints[i].type)
  1529.             zpoints[i].type = ypoints[i].type;
  1530.  
  1531.         }
  1532.         xicrvs = xicrvs->next;
  1533.         yicrvs = yicrvs->next;
  1534.         zicrvs = zicrvs->next;
  1535.         }
  1536.  
  1537.         /* Ok, fix up the title to include xp and yp plots. */
  1538.         if (((xp->title && xp->title[0] != '\0') ||
  1539.          (yp->title && yp->title[0] != '\0')) && zp->title) {
  1540.         tlen = (xp->title ? strlen(xp->title) : 0) +
  1541.             (yp->title ? strlen(yp->title) : 0) +
  1542.             (zp->title ? strlen(zp->title) : 0) + 5;
  1543.         new_title = gp_alloc((unsigned long) tlen, "string");
  1544.         new_title[0] = 0;
  1545.         if (xp->title && xp->title[0] != '\0') {
  1546.             strcat(new_title, xp->title);
  1547.             strcat(new_title, ", ");    /* + 2 */
  1548.         }
  1549.         if (yp->title && yp->title[0] != '\0') {
  1550.             strcat(new_title, yp->title);
  1551.             strcat(new_title, ", ");    /* + 2 */
  1552.         }
  1553.         strcat(new_title, zp->title);
  1554.         free(zp->title);
  1555.         zp->title = new_title;
  1556.         }
  1557.         /* add xp and yp to head of free list */
  1558.         assert(xp->next_sp == yp);
  1559.         yp->next_sp = free_list;
  1560.         free_list = xp;
  1561.  
  1562.         /* add zp to tail of new_list */
  1563.         *last_pointer = zp;
  1564.         last_pointer = &(zp->next_sp);
  1565.         xp = zp->next_sp;
  1566.     } else {        /* its a data plot */
  1567.         assert(*last_pointer == xp);    /* think this is true ! */
  1568.         last_pointer = &(xp->next_sp);
  1569.         xp = xp->next_sp;
  1570.     }
  1571.     }
  1572.  
  1573.     /* Ok, append free list and write first_plot */
  1574.     *last_pointer = free_list;
  1575.     first_3dplot = new_list;
  1576. }
  1577.